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Abstract 

We study a partially ionized hydrogen plasma by means of quantum molec- 
ular dynamics, which is based on wave packets. We introduce a new model 
which distinguishes between free and bound electrons. The free electrons are 
modelled as Gaussian wave packets with fixed width. For the bound states 
the ls-wave function of the hydrogen atom is assumed. In our simulations 
we obtain thermodynamic properties in the equilibrium such as the internal 
energy and the degree of ionization. The degree of ionization is in good agree- 
ment with theoretical predictions. The thermodynamic functions agree well 
with results from quantum statistics for 10000K ~ T ~ 40000K. 
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1 Introduction 



An understanding of dense plasmas is important to many areas of physics such as 
astrophysics or solid state physics. Because a complete analytic theory has not yet 
been advanced, several simulation techniques have been developed |TJ, El |3], |], ^ in 
order to determine properties of plasmas. In such systems, however, the fermionic 
character and the long-range Coulomb forces cause many difficulties. 

Quantum molecular dynamics (QMD) has been used to solve these problems 
|], [7j]. This method has been developed as an extension of classical molecular dy- 
namics ||. In QMD, wave packets are used as an approximation of the wave function 
of particles. Their dynamics can be derived from a time-dependent variational prin- 
ciple H |9|. It is possible to consider bound states, exchange effects and correlations 
in this method. 

The traditional QMD technique is based on Gaussian wave packets. The first 



simulations of this kind were carried out by Heller [10]. This technique was fur- 



ther developed in nuclear physics to study scattering processes of nuclei || |7j and 



recently has been successfully applied to plasma physics by Klakow at el. [||, [TT 
Both equilibrium and non-equilibrium properties of dense plasmas can be studied. 

In this paper we describe a new QMD model for partially ionized plasmas, which 
extends the previous work by incorporating free and bound states and transitions 
of the electrons. In plasmas, such transitions can be caused by radiation and by 
collisions. Only the latter transitions are important in the range of temperature and 
density that we are considering. The description of such processes is very compli- 
cated, especially because the microscopic treatment of the 3-particle recombination 
is still an open question. 

There have been recent attempts made to introduce transitions in QMD. In these 
it has been necessary to upgrade to Hamiltonian dynamics by use of additional el- 
ements. Ohnishi and Randrup have proposed to allow transitions on the basis of a 
stochastic dynamics |12|]. Tully has developed a branching concept in order to study 
electronic transitions at surfaces |TB[. In his approach the Hamiltonian dynamics 
was interrupted at certain points and a quantum mechanical transition was carried 
out. 

In our model we start with a similar idea. In order to describe ionization and 
recombination processes we implement two possible transitions: A free electron can 
ionize an atom in a collision, and an ion and an electron can recombine in a 3- 
particle collision of two free electrons and an ion. We use the well-known cross 



section for the ionization process fT3|, ITS]. In order to model the rather compli- 



cated recombination process we have developed a preliminary model. We study the 
partially ionized, weakly degenerate (n e A 3 <0.15) hydrogen plasma at a density of 
1.35xl0 22 cm~ 3 in the temperature range from 10000K to 200000K (see fig. [I]). This 
is a much lower density than those studied by Pierleoni et al. and Klakow et 



al. 1 11]. Simulations based on this model converge into an equilibrium of ionization 
and recombination with the degree of ionization in good agreement with theoretical 
predictions. 
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Figure 1: Characteristic lines of the hydrogen plasma in the temperature-density 
plane and data points of different investigations: • indicate our calculations, a are 
those of [16] and o are those of [11]. 

2 Basic assumptions 

2.1 Quantum molecular dynamics 

Quantum molecular dynamics is a computational technique to study quantum many 
body systems. It is based on an approximation for the single particle wave function. 
Such a trail state *f? qu (t) contains the essential degrees of freedom and the parameters 
q v (t) represent the coordinates in a generalized phase space. The time evolution of 
the trail wave function is specified by the dynamics of the time-dependent parameters 
q u (t). For a chosen parameterization of the trial state, the equations of motion for 
the parameters can be derived from a time-dependent variational principle || : 



dtlV 



.d ~ 
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dt 







(1) 



The equations of motion can be written in the general form 

{ -idH 



(2) 



in which the Hamilton function H is defined as the expectation value of the Hamil- 
tionian H, 
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is a skew-symmetric matrix JO, |I7|. In certain trial wave functions, the parameters 



can be grouped into pairs of canonical variables and the matrix assumes the canonical 
form 

I 1 



A 



(5) 



which leads to a Hamiltonian dynamics for the parameters. Gaussian wave packets 
(GWP) obey this type of dynamics. 



2.2 Free electrons 



GWP have a long tradition in QMD [|, |5|, |10|, p], fL8|, |7| . We use them as trial wave 
function to describe the free electrons in the plasma: 
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(6) 



The position r and the momentum p stand for the expectation values of the rele- 
vant quantum mechanical operators. The parameter 7 determines the width of the 
wave packet. Klakow et al. made the first plasma simulation with GWP. In 
their model 7 is a complex parameter of the dynamics and an additional degree of 
freedom of the trial wave function. We have shown that this ansatz leads at high 
temperatures to a heat capacity per electron greater than [17|. 

Because of this principle problem we fixed the width of the GWP in our simula- 
tions. A condition for 7 can be derived from the Kelbg potential |L9| : 
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which should approximately coincide with the interaction potential of two GWP [|TT 

e 2 
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Assuming V e e ) GWP-GWp\ n = V Kelbg we fitted 
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This ansatz leads to canonical equations of motion for the parameters r and p. 
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2.3 Bound electrons 



In our simulation we consider the ls-ground state of hydrogen explicitly. The best 
approximation for the ground state with GWP leads to an energy of 11.5 eV instead 
of 13.6 eV [ n|| . GWP with a dynamics derived from (|l|) have a continuous spectrum 



of excitation and a gap is missing. This results into unphysical excitations in many 
particle simulations and strongly influences the thermodynamic properties of the 
system |T7[. That is why we introduce a Is wave function (lsWF) to describe 
bound states, 
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where fj is the position of the core. This wave function is an eigenstate of the 
Schrodinger equation. 



2.4 Simulation of Hydrogen Plasma 

The hydrogen plasma consists of ions with positive charge and electrons with neg- 
ative charge. In order to incorporate quantum mechanical effects, the electrons in 
our model are treated as wave packets. The ions, which have much greater mass, are 
considered to be classical particles. We distinguish between free and bound states 
of the electrons by representing free electrons by GWP © and bound electrons by 
lsWF (fl0|). Our method is based on the Hamilton function H, which consists of the 
several parts, 

H = T i + T e + V ii + V ee + V ei . (11) 
The classical ions contribute kinetic and potential energy, 
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where the Ewald potential is used for periodic boundary conditions [T], |j. The 
following terms depend on the trial wave functions of the electrons, 
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where r is the distance between the particles. Different wave functions lead to 
different potentials for the electron-electron interaction. The interaction of two 
GWP is given by The interaction of a lsWF with either a GWP or a second 
lsWF is given by 
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An atom in the simulation consists of a proton and an electron in the ls-ground 
state. The interaction potentials of atoms are obtained by adding the contributions 
from both parts. 

The ansatz (^|) reveals the extra term 3/(2m7 2 ) in the kinetic energy of a GWP 
(pO"D , which has its origin in the shape of the GWP. Since T e ^ GWP is not consistent 
with the classical limit of free particles for high temperatures p 2 /2m we omit this 
term. 



3 Thermodynamic Model 

In this paragraph it will be shown how the QMD can be used to determine the degree 
of ionization z (the average charge number) of the plasma. In the usual approach 
one derives the degree of ionization from the free energy F, which is minimal at 
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constant volume and temperature. Alternatively it can be determined by using the 
condition that the internal energy U is minimal at constant entropy and volume. 




T.V,n 




(26) 



The free energy is not available in QMD simulations. That is why the second 
condition is used to find the equilibrium-degree of ionization. In order to find the 
minimum of the internal energy U a number of separate simulations with different 
degrees of ionization were carried out without including transitions between free and 
bound states. 

The entropy must be the same in all of these simulations. This requires a different 
temperature in every simulation. It can be obtained from entropy function S, which 
has to be known in our method. We use the entropy of the ideal hydrogen plasma, 
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iVfc is the number of particle in the volume V and g k is the degeneracy parameter 
due to the spin, where the index k runs over all particle types. If one describes a 



hydrogen plasma in the chemical picture pO one has to consider atoms (go=2), ions 
(<7i=l) and free electrons (g e =2). The degree of ionization is defined by z = N e /N, 
where N e =Ni and N = N + A^. 

Under the condition of fixed entropy the appropriate temperature can be de- 
rived by equation (p7|) for any degree of ionization. In the simulations we use a 
velocity-scaling procedure [f|, which forces the system towards the desired temper- 
ature. Then the internal energy is obtained by averaging over many phase space 
configurations of the system. In this way the internal energy can be determined 
point by point on an adiabatic and its minimum can be found. The exact value 
for the minimum is obtained by quadratic approximation (fig. |2|). The temperature 
belonging to this minimum has to be calculated from (^7|) since the temperature in 
all simulations is different. 



Results 



We study a partially ionized hydrogen plasma at a density of 1.35x10 cm in 
the temperature range from 10000K to 200000K. The plasma is weakly degenerate: 
n e A 3 < 0.15, where n e is the density of the free electrons derived from the ideal 
Saha equation [20|. We consider 32 electrons and 32 ions in our simulations. 

In figure ^, the degree of ionization from the simulation in the thermodynamical 
model is compared with results from the Pade approximation in the chemical picture 
(PACH) |f2TJ and from the ideal Saha equation. The PACH approach is based on 



Pade formulae for the thermodynamic functions and on the non-ideal Saha 
equation [20]. All methods show no ionization below 20000K, then the degree of 



ionization increases sharply. The values from our simulation and from the PACH 
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Figure 2: Internal energy u as function of the degree of ionization at constant 
entropy. 

are higher than that of the ideal plasma because the Coulomb interaction leads to 
a reduction of the effective ionization energy. At higher temperatures, the Coulomb 
interaction is less important, therefore the difference between the results from our 
simulation and the ideal Saha equation becomes smaller. The values from the PACH 
are still 6% higher because the Planck-Larkin partition function, which takes into 
account excited states is used in these formulas PP[ . 

The internal energy per particle from the simulation in the thermodynamic 
model u t h, from the PACH upach and from the ideal plasma are shown in figure 
|] and [5]. Under 20000K hydrogen behaves like an atomic gas, so u increases like 
3/2 kgT. All three theories coincide in this region. With increasing temperature u t h 
and upach rise more quickly than because of the higher degree of ionization, 
which leads to addtitional contributions to the kinetic energy. At high temperatures, 
negative contributions from the Coulomb interaction determine the corrections to 
Uid- So u t h and upach intersect with at an intermediate temperature. 

The agreement between the u t h and upach is fairly good in the region 10000K~ 
T~ 50000K. For higher temperatures, deviations start to occur. Above 75000K, uth 
and upach show the same behavior but differ by a constant. The PACH predicts a 
smaller energy although the degree of ionization is higher. One possible explanation 
is that there are no interactions of charged and neutral particles taken into account 
in these formulas. Further we have to mention that our approach neglects excited 
states, which are taken into account in the PACH approach |Tfl . 
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Figure 3: Degree of ionization versus temperature 



4 Dynamic Equilibrium of Ionization and Recom- 
bination 

There have been several attempts to model reactions by QMD. Tully developed 
a branching concept for transitions between different electronic states ||13|| . In this 
approach, the transition probabilities were derived from a Greens function technique. 
Hamiltonian dynamics is still used most of the simulation. Only in the case of a 
transition, dynamics are interrupted, new initial conditions for the particles are 
formulated and then the dynamics are continued. 

In our simulation we must deal with ionization and recombination processes. 
If an electron with an impact energy greater than 1 Ry collides with an atom, 
ionization can take place. In order to include such transitions in QMD, knowledge 
of the microscopic dynamics of such a process is necessary. The ionization of one 
hydrogen atom in a plasma was studied in p3 |. Here we investigate the much 



more complicated case where many atoms are present and include ionization as 
well as recombination. The ionization cross section for hydrogen is known various 
theoretical and experimental predictions WB, . We use the semi-empirical formula 



15, 23 
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Adapting the cross section to the wave packet model we imagine a sphere of 
radius ^ai(E)/n around every atom. If a free electron enters this sphere with an 
impact energy E greater the 1 Ry the ionization of the atom occurs. The electron 
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in the Is ground state is replaced by free GWP, which is place at the opposite point 
on the sphere. The momenta are chosen arbitrarily under conservation of energy. 

A microscopic description of the recombination process is more difficult. How- 
ever special information on recombination is available from various investigations. 
Besides the cross sections in different approximations which were derived from scat- 
tering theory, also the global recombination rate f3 is known from the theory of the 



rate coefficients fll4| , p24fl . A microscopic description of the recombination process 
in QMD has not been derived yet. As a preliminary model we propose the follow- 
ing mechanism: If two electrons are simultaneously in a sphere of radius Pr — I 
&b around an ion the recombination will take place. One electron is moved into 
the ground state and the other changes its momentum so that the total energy is 
conserved. 



Results 

We carried out simulations the the same density from 30000K (128 ions and 128 
electrons) to 200000K (32 ions and 32 electrons). The system converges into a 
dynamic equilibrium of ionization and recombination. The degree of ionization we 
get from this simple model is surprisingly good (fig.^J). It agrees well with the PACH 
approach in the considered temperature range. We do not claim that we give with 
our simple microscopic model a rigorous description of the recombination process 
but, rather, we deduce that our mechanism leads to a correct global rate (3 in the 
temperature range studied here. Correct in context means that the average number 
of recombination processes per time unit is in agreement with the number ionization 
processes so that the detailed balance is satisfied. 

In our model, we have found a dependence of [3 on pn and a setting of pn = lao 
turned out to be a good estimate. 

The internal energy u derived from this simulations leads to qualitatively correct 
results (fig.fD but quantitatively the results only agree well up to temperature of 
40000K (fig.[|). Beyond of that our results are systematically high in comparison to 
the PACH, which may be caused by an insufficient description of the electron-ion 
interaction, because our results were generally too high in comparison to the PACH 
model. Evidently simulations above 40000K require a more precise description of 
the interaction between free electrons and ions than that given in our QMD model. 



5 Conclusion 

We expect that the model presented in this paper will prove to be a promising 
approach for the description of plasmas in the partially ionized region. We have 
substantially extended the traditional QMD model with the consideration of free 
and bound states. We have shown how thermodynamical functions can be derived 
from a model with free and bound states but still without transitions. Furthermore 
we presented preliminary model for the dynamic description of ionization and re- 
combination in plasmas. As first test we calculated the degree of ionization and the 
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internal energy of the hydrogen plasma. 

The model can be improved in several ways. First of all a detailed quantum me- 
chanical description of the 3-body recombination process would be useful in future 
studies. Furthermore it would be interesting to study the formation of molecules 
One must also take into account the fermionic character of the electron wave 



function J7|, |TT|| . Our results offer hope that it soon will be possible to describe 
hydrogen in the whole temperature range from the molecular gas to the completely 
ionized plasma. 
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